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Abstract 

Magnetic stimulation is a standard tool in brain research and many 
fields of neurology, as well as psychiatry. From a physical perspective, one 
key aspect of this method is the inefficiency of available setups. Whereas 
the spatial field properties have been studied rather intensively with coil 
designs, the dynamics have been neglected almost completely for a long 
time. Instead, the devices and their technology defined the waveform. 

Here, an analysis of the waveform space is performed. Based on these 
data, an appropriate optimisation approach is outlined which makes use 
of a modern nonlinear axon description of a mammalian motor nerve. The 
approach is based on a hybrid global-local method; different coordinate 
systems for describing the continuous waveforms in a limited parameter 
space are defined for sufficient stability. The results of the numeric setup 
suggest that there is plenty of room for waveforms with higher efficiency 
than the traditional shapes. One class of such pulses is analysed further. 
Although the voltage profile of these waveforms is almost rectangular, 
the current shape presents distinct characteristics, such as a first phase 
which precedes the main pulse and decreases the losses. The single rep- 
resentatives, which differ in their maximum voltage shape, are linked by 
a nonlinear transformation. The main phase, however, seems to scale in 
time only. 

Keywords: Magnetic stimulation, TMS, inductive stimulation, waveform op- 
timization, neuron model, loss power. 
PACS: 87.19.1b, 87.19.1d, 02.60.Pn, 87.50.fc 

1 Introduction 

Magnetic stimulation is a standard tool for noninvasive activation of neurons, 
especially in the brain (1). The underlying physical principle of this method is 
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induction, which is independent from any physical contact. Due to the magneti- 
cally passive behaviour of biological tissue, the magnetic field of the stimulation 
coil can pervade the body and especially poorly-conducting structures, such as 
bones. For a supra-threshold stimulation of brain neurons using electric cur- 
rents, in contrast, the skull acts as a barrier. The consequently higher currents, 
however, lead to distress which is usually not tolerated by patients without 
anaesthesia (2). Moreover, the current paths are no longer obvious, whereas 
inductive stimulation provides a relatively focused method for local stimulation 
of small neuron populations. 

However, magnetic stimulation is extremely inefficient. Less than one per- 
cent of the electrical energy is transferred to the target. Although a substantial 
fraction of the field energy can be recuperated and used in the next pulse, the 
high currents in the range of kiloamperes during operation lead to significant 
ohmic losses in cables, connectors, and the stimulation coil itself. Driving such 
high charge flows necessitates up to 4000 V. The application of high voltage is 
a serious issue for patient safety; the high losses heat the coil and limit typical 
durations of repetitive stimulation sessions to several minutes. Mobile devices 
with battery-powered pulse sources were proposed, but are not reasonable in 
the context of repetitive protocols. 

Stimulation of neurons comprises two important domains. The current of the 
stimulation coil induces an electric field in the tissue, which in turn leads to drift 
currents. The role of induction with respect to efficiency was already discussed 
before (3, 4, 5). Two aspects lead to the small energy transfer in induction. The 
low conductance of biological tissue, which is more than six orders of magnitude 
lower than in copper, counteracts higher eddy currents. The second issue is the 
coupling between the primary windings — the stimulation coil here — and the 
tissue or a distinct region thereof, which corresponds to the secondary side of 
this transformer. Improving the coupling is principally a geometric task, which 
may incorporate magnetic materials (6) . This spatial degree of freedom can be 
covered with appropriate coil designs (7, 8, 9, 10, 11, 12). 

The key to the temporal perspective is the second stage, namely the phys- 
iology of the neurons. Whereas the induction process is approximately linear, 
the neuron dynamics introduce a complex nonlinear system that may not be 
neglected in this context. In the history of magnetic stimulation devices, ef- 
ficiency has played a secondary role for a long time. First systems primarily 
aimed at reaching the stimulation threshold noninvasively with the available 
technological means — within more than a hundred years since J. C. Maxwell, 
many approaches have failed until this goal was finally reached (13, 14, 15, 16). 

In the first devices, the entire energy was lost during a pulse. This so-called 
monophasic topology had been the standard for many years, until the biphasic 
oscillator circuitry introduced some improvements due to omitting the inten- 
tional damping. These were driven by purely technical considerations and al- 
lowed feeding back a notable fraction of the magnetic field energy into the pulse 
capacitor. The 'needs' of the neurons were not taken into consideration also in 
the later device generations. Originally, it was thought that the falling phase 
of such an oscillation might counteract the activation of the leaky neuron mcm- 



2 



brane (5, 17). However, the system turned out to be not only rather efficient, but 
even more effective for the same voltage amplitudes (11, 18, 19, 20, 21, 22, 23). 

The question of dynamics is already an older issue in electrical stimulation. 
A lot of parameter studies were performed with certain predefined waveforms 
(24, 25, 26, 27, 28, 29, 30, 31, 32). An analytic optimisation study, which led 
to the so-called rising exponential pulse, was reported already 1946 based on 
a variational approach (33). A lot of work in this field has been done with 
Lapicque's well-known linear model from 1907 (34). However, shortly after 
Offner's optimisation approach, it was shown that neuron dynamics are based on 
the interplay of a vast number of single components which are highly nonlinear 
(35). Despite that, his work with linear dynamics was reproduced over sixty 
years later (36), but the results are neither compatible with more sophisticated 
neuron models (37) nor with experiments (27). 

Nonlinear models, in contrast, can no longer be inverted easily, and chaotic 
behaviour of the corresponding differential equations render also the numerical 
handling problematic. A first general, (nearly) unbiased optimisation of the 
waveform for electrical stimulation was done recently (38). The authors min- 
imise the ohmic losses of monophasic electrical pulses. The optima reported 
there seem to be far more likely Gaussian than rising exponential. 

For magnetic stimulation, there are only some rare parameter studies avail- 
able in the literature related to different waveforms using experiments (39, 23, 
40, 21) and numerical models (41, 42, 43, 44). Thus, even the knowledge about 
the different characteristics of existing waveforms is limited at the moment. 

The question of optimality has not been studied so far. Instead of discussing 
again any technological changes of the device and their effects on neurons, this 
text wants to change the perspective and tries to address the requirements of 
a neuron instead, without any unnecessary restrictions, using a more realistic 
model. 

2 Setup 

2.1 Objectives/ Aims 

The central aim of this approach is a general optimisation without unnecessary 
constraints. Thus, no limitations of existing devices are taken into account. 
Only physical principles are of importance here. 

For the objective function, several aspects of magnetic stimulation can be 
used. Here, a key problem of the technology is minimised, namely the energy 
losses. The heating of the stimulation coil limits the maximum duration of a 
session, as well as the repetition rate (45). Furthermore, the losses have a direct 
impact on the power consumption. The fact that portable repetitive devices are 
still missing is a consequence of that. 

If no specific losses of a specific technology, such as the forward bias of p-n 
junctions in a concrete circuitry, are regarded, the main cause of power drain 
that can hardly be avoided is Joule heating relating to the inner resistance of 
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all components. For a coil current i(t), the losses of a pulse regardless of the 
exact value of the inner resistance Ri ^ — are J R+ Rii 2 dt cx J R+ i 2 dt =: W. 
This integral forms the objective functional W for the further analysis. 

In addition to the objective, several constraints are of importance — even 
without the limitations of a certain technology. Firstly, the coil current has to 
elicit a neuron response. A second constraint results from the interplay of the 
objective and the response condition. As will become apparent later — and is 
known to many device designers — the losses at the excitation threshold fall with 
shorter pulse durations for classical waveforms. At the same time, the required 
device voltage increases (42). In a consequent optimisation approach which 
does not regard that, the voltage diverges consequently. This degree of freedom 
can be determined by several means from a physical perspective. A simple 
voltage limit v max which may not be exceeded by the coil voltage \v c \ < v max in 
positive as well as in negative polarity was chosen here. Although this is often 
motivated by the limitations of semiconductor devices — regardless of the exact 
technology — it might be most adequate also with respect to patient safety and 
insulation. For being independent from the exact inner resistance, the voltage 
v c relates to the voltage of the inductance of the coil only. The current profile, 
in contrast, is independent from the losses. 

2.2 Framework 

The optimisation framework has to describe the temporal behaviour, including 
all stages from the coil through to the response. Figure 1 depicts the structure 
of the single modules. The setup is outlined shortly as far as this information 
supports the understanding of the results. 

The coil current induces an electric field. From a dynamic perspective, in- 
duction is a differentiating process. It is assumed that the tissue around the 
axons does not contain too many very different material types so that the filter- 
ing and distorting influence is rather low, and the first-order derivative becomes 
dominating. The latter provides the excitation function for the axon model. 
Further terms for describing distortion effects are optional, but set to here. 
The neuron model, in turn, hands back the first constraint. 

The second constraint and the objective function relate to the device as 
mentioned. Their evaluation is performed as described dependent on the voltage 
v c and the current i. 

2.3 Parametrisation/Problem Description 

The objective and all constraints were formulated as functionals. For numerical 
optimisation, however, this abstract construction has to be translated into a 
set of numbers. Thus, an appropriate definition of a parametrisation is an 
essential element, which moreover has a notable influence on the stability of 
the optimisation system. This mapping from the finite parameter space to a 
subspace of the Hilbert space of all waveforms allows the stable evaluation of 
all needed operators, such as derivatives, in an analytical way. Evading this 
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Figure 1: Structure of the optimisation setup. The minimisation system itself 
has a hybrid structure. A global framework commands many local workers. 
The algorithms minimise the weight W in the space of valid constraints C, but 
has only access to a finite number of abstract parameters p. The latter are 
converted to functions by ^} as described. The initialisation is performed with 
noise, as well as results from earlier runs. The axon model is incorporated into 
the constraints as a nonlinear element (AfC). 

issue by taking a large number of sample points which are handed over to an 
optimisation algorithm (38) might end up in a lengthy brute-force approach and 
instabilities. 

The outlined conditions of magnetic stimulation are rather complex. Espe- 
cially induction with its differentiating behaviour intensifies the problem. Lots 
of different waveforms are able to stimulate a nerve — their common features, 
however, are hardly recognizable. The part of the surface which is formed by 
the objective and furthermore fulfils the first constraint shows a strong tendency 
to a very rugged 'crater landscape'. 

The selection of an appropriate coordinate system for the description of the 
waveforms can remarkably improve the solution process, but the problem is non- 
linear and has no natural (discrete) parametrisation. Three simple coordinate 
systems were used simultaneously here. 

Firstly, the waveform is described by cubic spline curves; their parameters 
p € K™ act as degrees of freedom for the optimisation algorithm. Accordingly, 
this approach is similar to the core principles of the finite element method in 
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numerics. The predefined solution forms a subspace of the class of all func- 
tions with continuous second derivative C 2 (T) with T being the limited time 
domain of the reduced problem. However, in view of inductive stimulation, this 
coordinate system still has lots of unfavourable local minima in the objective 
function. 

A parameter description in the frequency space was taken as an alternative. 
A complex-valued parameter vector p' e <C n / 2 can be taken without further 
constraints. A discrete Hilbert transformation provides the parameters for the 
single terms of a (complex-valued) Fourier series. Accordingly, the latter is even 
a smooth junction and a member of C°° . 

Because the here applied optimisation algorithms work with simple floating- 
point numbers only, there are two ways for generating the complex-valued pa- 
rameter vector p' from a real- valued vector p € K n according to standard math- 
ematics. On the one hand, the vector can be split into a real and an imaginary 
part. On the other hand, a partition into the absolute value and the phase is 
possible. 1 

All of the three proposed methods of parametrisation have at least two con- 
tinuous derivatives. All of them allow a dynamic change of the coordinate 
system — this comprises a conversion from one parametrisation to another, as 
well as increasing or reducing the degrees of freedom in the same type of de- 
scription. All of them were incorporated here. 

This may look like indecisiveness, but it was done intentionally in order 
to address the main problem of this endeavour, namely instability due to a 
vast number of local minima. A local minimum in one space might not be so 
pronounced in another. Lots of observed local minima showed artefacts typical 
for one type of description, such as Gibbs phenomena or ringing, which are not 
stable for the optimiser with another parametrisation. A typical example for 
a whole class of waveforms are the harmonic functions. A sinusoidal current 
stimulates as does its negative, i.e. mirrored version. In the time-domain, the 
shortest transition in-between passes 0, which does not stimulate at all and 
separates two minima. In certain frequency-domain descriptions, this is a less- 
problematic continuous shift of the phase parameters without such barriers. 
Furthermore, a dynamic change of the number of dimensions was used for an 
adaptive control over stability. For a refinement, the degrees of freedom were 
increased. A diverging run can be re-stabilised by reducing complexity. 

The parameters describe the current shape. Although a parametrisation of 
the electric field could avoid the differentiation of the function — integration is 
required instead — the results did not improve in that case. The duration of the 
support of the waveforms was reduced to three milliseconds during the piloting 
phase of the setup as no local minimum was observed to require a longer time 
base. 

In case of the second method, the phase values can be increased by a constant factor 
in order to match their magnitude with the absolute values. This is an important issue for 
simplex-based optimisation methods as were used here. Those estimate the deviation usually 
from the next minimum on the basis of the area of a simplex. In the case of COBYLA, the 
goodness of the simplex additionally decides over the type of the next step. 
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2.4 Optimisation Algorithm 



It seems a noteworthy hint in this context that not the objective, which is just 
squared, is the most challenging functional here, but the first constraint. 

An optimisation method for multivariate minimization of the parametrised 
objective was implemented. Due to the high number of local minima, a global 
method is advisable. For a quick convergence, this was combined with a local 
algorithm in a hybrid approach. 

Local optimisation workers with their much faster convergence run into the 
next local minimum. A global framework in turn combines the information 
about local optima. 

For the global algorithm, a particle swarm was selected (46), which seems to 
be more systematic in this context than defining genes for certain combinatory 
operations in a genetic algorithm, for instance. In addition to that, a particle- 
swarm framework with lots of local workers very easily allows a concurrent 
design that meets the requirements of high-performance computing. 

Two algorithms — a simplex descendent (constrained optimisation by linear 
approximation, COBYLA, (47)) and an interior-point method (48) — act alter- 
natively as relatively stable local optimisers. The type of the coordinate system 
of a certain worker is fixed. The number of degrees of freedom applies to all 
workers and is controlled by the global framework. The number is increased 
by a predefined step if convergence is achieved and the result outperforms the 
latter; otherwise, it decreases. 

The global method hands over the parameters to the local workers, which are 
supposed to run into a nearby valid local minimum which fulfils all constraints. 
Their results are regarded for deciding on the particles' as well as the total best 
in each step. In order to avoid oscillations in the attraction field of local minima, 
those are not assigned to the current position of a particle. 

The update rules of the particle swarm assign the following to the parameters 
of the j-th particle in the (i + l)-th step 



with the inertia w, the gravity parameters ci, ci of the local and the global 
best, as well as the modulation variables r^^ji € (0, 1). The latter are chosen 
randomly in every step. The global best is denoted by b ff , whereas each particle 
has its own local best hj. For the constants, several alternatives were tested. 
The following values performed appropriately: ui = 0.6, c\ = 1.7, ci = 1.7, 
alternatively even with repellent behaviour according tow = —0.35, c\ = —0.05, 
c 2 = 5. The number of particles was up to fifty. 

Initialization of the single particles is done with random numbers. In several 
runs, a fraction of the workers was initialized with the known classical wave- 
forms. 
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2.5 Neuron Model 



The decision if a pulse elicits an action potential is performed by a nonlinear 
model of a human motor axon. The local description (see the appendix) , which 
is an adaptation from the literature, includes fast sodium channels, persistent 
sodium channels, one type of potassium channels and passive leakage assembled 
on the basis of (49, 50, 51, 52, 53). Its consistency with stimulation experiments 
which address especially the nonlinear behaviour of nerves was studied in (54). 
A response that exceeds +10 mV is seen as successful. 

The whole implementation is done in C and optimised for maximum speed. 
For the solution of the differential equation, a standard second-order Runge- 
Kutta method was used. The allowed maximum time step for the final analysis 
was set to 500 ns. The whole optimisation framework followed a strict parallel 
design. Every call of the differential-equation solver is executed in an own 
thread. The computation was performed on three Xeon servers with eight cores 
in each for the pre-evaluation and tests, as well as on a high-performance system 
of the Leibniz Supercomputing Centre of the Bavarian Academy of Sciences and 
Humanities. The results of this paper are based on more than 50 000 CPU hours. 

3 Pre-Evaluation 

The study of different waveforms from a physical perspective has been neglected 
so far in inductive stimulation. Even basic relationships are not generally known, 
but are mostly specific experience of academic and commercial device designers. 
Compared to electrical stimulation, the circumstances are more complex because 
of induction. Therefore, they are also less obvious. 

For that reason, a pre-evaluation of solutions was performed. Additionally, 
this acted also as a test of the stability of the algorithms at the same time. The 
results are presented in figure 2. This graph depicts the losses in dependence 
of the maximum required voltage at the threshold. All values are relative be- 
cause the threshold is nevertheless an individual figure. However, the voltage is 
roughly in the range of real devices if the unit volt is added. Exact calibration 
can be conducted easily based on the data. 

Every cross in the graph represents a local optimum. In addition to them, 
lines for the standard biphasic pulse (red line) and rectangular voltage pulses 
(cyan, dashed) were added from a parameter study for a frequency range from 
500 Hz to 50 kHz. Moreover, the waveform of a commercially-available monopha- 
sic device (Magstim 200, Magstim Ltd., Whitland, Wales) was added (star). For 
all of them, that current direction was taken into account which presented the 
lowest threshold. 

Most interesting in this graph is that the monophasic waveform is not notably 
less efficient than a biphasic pulse in the model. It seems to be the way of 
generation with the intentional damping only that wastes the pulse energy. This 
is remarkable because the current tail is rather long so that it makes a substantial 
contribution to the loss integral. This might call for better topologies since 
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Figure 2: Local minima from the pre-evaluation plotted in a plane which is 
spanned by the peak voltage and the pulse losses at the excitation threshold. 
Additionally, biphasic and rectangular (with respect to the voltage) pulses were 
depicted for a wide range of carrier frequencies. The insets outline the current 
shapes some examples. The small letters refer to solutions which are shown in 
figure 4. 



monophasic pulses have several distinct features, such as the higher sensitivity 
to the coil orientation, which is favourable in a lot of applications (19, 55). 

Nearly all waveform types seem to have the general tendency that the losses 
can be reduced if a higher voltage is provided in the device. This was al- 
ready mentioned in the text above. For biphasic and rectangular stimuli, the 
transformation is just a simple dilation or compression. Although this effect 
is sometimes used for reducing the heating in high-power applications, e.g. in 
rehabilitation (56), there seem to be pulse shapes that have a much steeper 
slope in this loss-voltage plane, i.e. the reduction of the losses for an increased 
voltage is notably higher than for biphasic pulses. 

Rectangular pulses might be more efficient than the classical waveforms. 
This assumption was already stated before (57, 39). However, they are still far 
from being optimal — especially for the voltage range of typical devices. 
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Figure 3: Relative current and voltage profiles for one exemplary minimum. 
The absolute values on the axis depend on many parameters, such as coil design, 
inductance, distance, etc. The full-width at half maximum of the main phase of 
the current pulse amounts to 23 us (see figure 2(c)). The maximum derivative 
within the first phase is by about a factor of 161 lower than the slopes of the 
second. Despite the low amplitude, the long duration causes that the area under 
the curve for the first current phase is higher than that of the second (by 15 %). 



The lower edge of the point cloud is smooth and seems to represent the same 
class of stimuli. For very high voltages, the current waveforms degenerate to 
Dirac-like pulses. 



4 Optimisation Results 

4.1 Basic Features 

For several different voltages in the range of classical devices, a more exact 
analysis was performed including a computationally- intensive refinement. The 
number of degrees of freedom was increased up to one thousand during these 
runs. Most of the computational resources were spent on this investigation. 
Some exemplary results are depicted in figures 3 and 4. 

The best found solutions are very similar. The current always exhibits three 
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relatively distinct phases. A slow first phase precedes a short main pulse in op- 
posite direction; a separate third phase becomes visible after very high iteration 
numbers of the optimisation algorithms. 

The shapes of the slopes within the single phases seem to have only minor 
influence on the losses. However, they define the maximum voltage level. The 
fine-analysis reveals an almost rectangular voltage shape of the second phase. 
Accordingly, the current is approximately triangular there. The other parts of 
the pulse seem to disappear in the voltage profile. 

The voltage, however, cannot be the benchmark because the pre-evaluation 
clearly shows that rectangular pulses are definitely suboptimal in this model. 
It has merely indirect influence. The current, in contrast, shows rather subtle 
characteristics. 

The square-shaped voltage during the main phase might be the consequence 
of the second constraint. The system uses the maximum allowed voltage level as 
long as possible. This explains the positive wing of the voltage. The negative, 
however, reduces the current as quickly as possible after the peak. Although this 
might counteract excitation, a lower current level — which has squared influence 
in the losses — might be more important due to the objective. This occurs only 
until a certain low level is reached from which the third phase heads for in a 
slower, exponential decay. 

This reasoning can also be used for making the first phase plausible. The 
latter biases the onset of the second phase with minimum dynamics. Start- 
ing from this shifted baseline, a rather long rising slope can emerge, without 
reaching extremely high loss powers. Also this first phase causes losses, but the 
amplitude is relatively low. Since the heating depends on the squared current, 
this investment seems to pay off. Prolonging the slope to the left, thus to nega- 
tive values, leads to much lower costs in the objective than increasing the peak 
current. For the waveform of figure 3, the losses would be approximately 11 % 
higher without the first phase; for lower voltages this difference is much higher. 2 

Nevertheless, it should be noted that such explanation is always an over- 
simplification. Ascribing certain functions to different phases of a pulse can 
be misleading in connection with nonlinear dynamics, although older literature 
might motivate such a step. 

4.2 Common Ground of the Corresponding Class of Wave- 
forms 

The single curves in figure 4 suggest that these waveforms can be treated as 
representatives of the same class. In contrast to biphasic or rectangular pulses, 
the common ground of the members is not a linear transformation, such as 
scaling. The behaviour of the first phase — which, for instance, even decreases if 
the main peak becomes stronger — speaks out against such a simplification. 

2 In the voltage range of the incorporated commercial monophasic device, the losses of the 
rectangular pulse seem to be even 25 % higher than in the corresponding local minimum with 
the same voltage limits (see figure 6). 
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Figure 4: Local pulse optima for different threshold-voltage constraints. The 
comparison reveals the different roles of the single parts. Whereas the length of 
the second phase roughly scales inversely to the voltage constraint, the duration 
of the first part stays constant. However, the amplitude of the first phase relative 
to the second varies. The exact data in the coordinate system of figure 2 are 
as follows: 334 as the corresponding loss value for a relative voltage of 990 (a), 
272 for a relative voltage of 1960 (b), as well as a loss value of 214 for a voltage 
of 4800 (c). 



If the second phase is isolated, however, a simple time scaling can be assumed 
in the first approximation. The pulse duration decides on the voltage level and 
influences the heating: For shorter triangles, the amplitude has to be increased, 
but the losses fall. The dynamics of the first phase, however, do not change at 
all for the studied range, but only the amplitude is adjusted for the different 
voltage constraints. 

The exact behaviour of such details was analysed in a parameter study. The 
amplitude of the first phase, i.e. the value of the most negative point in the 
current profile, was rendered changeable in dependence of the main-phase peak 
amplitude. As a second degree of freedom, the duration of the rising slope of 
the main phase was varied; the falling slope changes proportionally to the latter. 
The optimal amplitude of the first phase was evaluated for a number of different 
main-pulse durations. The results are presented in figure 5. 

The curves depict the amplitude relation as well as the contribution of the 
first phase to the objective, thus the value of the integral of the squared current 
over the first phase in relation to the total weight. For very short pulses, the first 
part of the waveform vanishes almost and the pulse becomes nearly triangular 
in the current profile, as well as rectangular with regard to its voltage shape. 
Also the investment in the prephase seems to be less rewarding. For the longest 
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Figure 5: Role of the first phase of the current for different representatives of the 
whole class of local minima. If the pulse duration is increased, the amplitude 
of the first phase rises relatively to the main pulse and leads to a nonlinear 
transformation rule between the pulses. The contribution of the first phase to 
ohmic losses is also not constant. The x-axis refers to the length of the rising 
slope. This is equivalent with the duration of the first voltage peak. 

depicted main pulse durations, however, the current at the apex of the first 
phase is nearly half of the main amplitude. 

Figure 6 shows the effect for different pulse durations in relation to sym- 
metric triangular current profiles, i.e. rectangular voltage pulses, with the same 
duration. With the increasing amplitude of the first phase in pulses with a 
longer second part, the advantage over the simpler competitor increases. For a 
duration of the rising triangular slope of 200 us, the local minimum amounts to 
62 % of the losses for the corresponding triangular current shape only. Again, 
for very short pulses, the difference in weight between the local optima and the 
triangular current pulses vanishes because also the differences in respect of the 
shape fade away. 

For the question of plausibility, a short look at corresponding results from 
another model might be worthwile. Taking the original Hodgkin-Huxley equa- 
tions and modifying especially their dynamics to mammalian body temperatures 
once was a common simple and moreover computationally quick approach for 
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Figure 6: Influence of the pulse duration on the advantage of the local minima 
compared to simple triangular current pulses. The y-axis depicts the relation of 
the ohmic losses of the corresponding local optimum to the heating of the cor- 
responding triangular current profile (hence rectangular voltage) with the same 
main pulse duration; both are at the excitation threshold. Even for unusually- 
short pulses with some tens of microseconds, the gain is still more than 10%. 

studying pulse dynamics, e.g. of the acoustic nerves (58). Although it is not 
appropriate for the analysis of different waveforms due to the notably different 
dynamics, the best found result from figure 7 also shows the known key phases. 

5 Conclusions 

In this text, a minimisation approach was formulated for inductive stimulation 
of an axon model; the key objective was the losses. Their minimum level is 
given by an ohmic linear voltage-current relationship, which leads to a squared 
dependence on the current for the heating. 

The approach was not limited to a specific circuitry, but the constraints were 
reduced to the necessary conditions only. The parametrised waveform space 
covers a very wide range of functions and is — from a practical perspective — 
almost unconfined. 

Although many of the worse local minima that emerged during the compu- 
tations presented rather curious pulse shapes, the best results found here seem 
very consistent. The single parts of the waveforms can be made plausible as 
shown, although this should not be misinterpreted as an explanation, which is 
a critical task in the context of nonlinear dynamics. 

All waveforms were reported with their current profile. This was done in or- 
der to emphasise the different roles of voltage and current for the predefined ob- 
jective. In the past, the voltage profile was usually foregrounded. Also in typical 
devices, the voltage level is controlled rather than the current (11, 12, 4, 59, 57). 
Concerning the losses, however, the current is the key figure. For the voltage 
profile, in contrast, the reference is not even well-defined. As the currents are 
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Figure 7: Local minimum from another axon model (Motz-Rattay) . The es- 
sential elements from the more realistic model above appear also here. 

relatively strong, the voltage drops at connectors, cables, etc. are high enough 
that the differences of the voltage shapes due to moving the monitoring points 
is notable. 

The potential of waveforms — not only for optimisation, but also in the con- 
text of additional degrees of freedom — seems to have only secondary prior- 
ity for commercial device manufacturers. The authors are not aware of any 
commercially-available system that would be able to provide the here derived 
shapes. This pulse generation has to be necessarily with recuperation of the 
field energy so as not to counteract the objective. Recently, a sophisticated 
technology was proposed by an academic team in order to provide more flex- 
ibility (57, 60). Although the described embodiment, which was designed for 
relatively low negative voltages, is not fully compatible with the aims given 
here, the underlying principle of that system could provide similar waveforms, 
at least. 

An alternative technology that would allow even a rather exact reproduction 
of such pulses was developed most recently (54) and will be discussed separately. 

From the perspective of efficiency, another tradition has to be reviewed. 
In classical devices, the amplitude is used for controlling the strength. Due 
to the 'monochrome' oscillator design, an alternative was not possible before 
cTMS (57). In the light of the falling losses for all classical pulse shapes, it 
could be advantageous in several high-power applications (56) to use the highest 
available voltage level of a device and modulate the strength with the pulse width 
instead. For sessions where the dynamics are important for certain effects, such 
as selective stimulation or diagnostics, the established amplitude control can be 
applied. 

Although the results were obtained in a systematic way and based on a 
plausible explanation, statements derived from any type of model should in 
general not be taken as irrevocable facts, but are rather a suggestion for ex- 
periments. The outcome was not biased by an anticipating initialisation with 
current or voltage shapes that resemble the final outcome, but the waveforms 
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evolved themselves. Many artefacts and worse local minima show much different 
curves. 

Despite that, the results are only (the best found) local minima. Moreover, 
they are relatively stable and re-emerged also in totally independent runs with 
random initialisation of the parameters. Indeed, these solutions can be global 
or not far from the global solution, but they do not have to. 

Furthermore, the found waveforms are undeniable local minima of the prob- 
lem formulation, which is based on a nonlinear axon model. Despite calibration 
and tests, each model shows limitations, shortcomings, and does only mimic 
reality. Especially in the absolute task of minimization, a physiologic optimum 
could look slightly different, and even be an individual quality of a nerve. Thus, 
this is also a request to experimentalists to further study this issue. 
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Appendix 



The state vector <j> of the nerve model is defined as follows. The first element 
denotes the membrane potential; the remaining items describe the channel dy- 
namics, which are represented by first-order systems with non-linear time con- 
stants: 
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c m = 2.0uF/cm 2 

g Na = 290mS/cm 2 
gNa, P = 25mS/cm 2 
g K — 80 mS/ cm 2 
gi = 7mS/cm 2 
E Na = 50.0 mV 
E K = -90.0 mV 
Ei = -90.0 mV 
T = 310K 

S rn = 1.86 (mVms)- 1 

e m = 21.4 mV 

C m = 10.3 mV 

T] m = 0.086 (mVms)- 1 

9 m = 25.7 mV 

= 9.16 mV 
S h = 0.062 (mVms)- 1 
e h = 114.0 mV 

= 11.0 mV 
rjh = 2.3 ms -1 
6 h = 31.8 mV 
L h = 13.4 mV 



S p = 0.01 (mVms)- 1 

t p = 27 mV 

C P = 10.2 mV 

r7 P = 2.5-10- 4 (mVms)- 1 

6 p = 34 mV 

l p = 10 mV 

S s = 0.3 ms -1 

e s = 53 mV 

Cs = -5mV 

rj s = 0.03 ms -1 

6 S = 90 mV 

L s = — 1 mV 

T r . act = 293K 
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for (pi = -e m 
else 

for cj)i = -6 m 
else 

for 0i = -e h 
else 

for 0! = -e p 
else 

for 0i = — 9 p 
else 
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